%--------------------------------------------------------------------------
%   小波分解重构工程实现
%   https://www.cnblogs.com/lianjiehere/p/4239342.html
%--------------------------------------------------------------------------
clear;clc;

%--------------------------------------------------------------------------
%   创建信号
%--------------------------------------------------------------------------
ns = sin(0:0.01:10-0.01);
ns = awgn(ns,30);

%--------------------------------------------------------------------------
%   信号进行尺度分解和小波分解
%--------------------------------------------------------------------------
%   matlab 自动分解
%--------------------------------------------------------------------------
[c,l]=wavedec(ns,4,'db2');                                                  %db1小波4层分解 分解矩阵+各层长度

%--------------------------------------------------------------------------
%   手动分解
%   获取db1小波系数
%   l_d     尺度函数分解系数
%   h_d     小波函数分解系数
%   l_r     尺度函数重构系数
%   h_r     小波函数重构系数
%--------------------------------------------------------------------------
[l_d,h_d,l_r,h_r]=wfilters('db7');
tempm=conv(ns,l_d);                                                         %尺度卷积
tempn=conv(ns,h_d);                                                         %小波卷积

filter_N = length(l_d);

m=tempm(filter_N:2:end);                                                    %掐头去尾 隔点抽取
n=tempn(filter_N:2:end);                                                    %掐头去尾 隔点抽取

%--------------------------------------------------------------------------
%   增采样
%--------------------------------------------------------------------------
m_up = dyadup(m);                                                            %增采样 插0
n_up = dyadup(n);                                                            %增采样 插0

m_up = m_up(1:1000);
n_up = n_up(1:1000);

%--------------------------------------------------------------------------
%   重构
%--------------------------------------------------------------------------
m_back=conv(m_up,l_r);
n_back=conv(n_up,h_r);

ns_back=m_back+n_back;
ns_back = ns_back(2:end-filter_N+2);                                        %去滤波器延迟 这里没搞懂为啥-2

figure(1);
plot(ns,'LineWidth',3);
hold on;
plot(ns_back,'-.','LineWidth',2);
plot(ns-ns_back);
hold off;title('对比')
